Magnetoelectric polarizability and axion electrodynamics in crystalline insulators 
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The orbital motion of electrons in a three-dimensional solid can generate a pseudoscalar magneto- 
electric coupling 8, a fact we derive for the single-particle case using a recent theory of polarization 
in weakly inhomogeneous materials. This polarizability is the same parameter that appears in 
the "axion electrodynamics" Lagrangian ACem = (6e 2 /2irh)E ■ B, which is known to describe 
the unusual magnetoelectric properties of the three-dimensional topological insulator (6 = tt). We 
compute 9 for a simple model that accesses the topological insulator and discuss its connection to 
the surface Hall conductivity. The orbital magnetoelectric polarizability can be generalized to the 
many-particle wavefunction and defines the 3D topological insulator, like the IQHE, in terms of a 
topological ground-state response function. 
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Magnetoelectric couplings in solids have recently been 
the subject of intense experimental and theoretical in- 
vestigations [3, 0, [H . A quantity of central importance is 
the linear magnetoelectric polarizability atj defined via 
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where E and B are electric and magnetic fields, P and M 
are the polarization and magnetization, and the equality 
can be obtained from commuting derivatives of an ap- 
propriate free energy. In general the tensor a has nine 
independent components, and can be decomposed as 



(2) 



where the first term is traceless and the second term, 
written here in terms of the dimensionlcss parameter 9, 
is the pseudoscalar part of the coupling. Here we focus 
on magnetoelectric coupling resulting from the orbital 
(frozen-lattice) magnetization and polarization, which we 
label the orbital magnetoelectic polarizability (OMP). 

In field theory, the pseudoscalar OMP coupling is said 
to generate "axion electrodynamics" [H , and corresponds 
to a Lagrangian of the form (c = 1) 
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An essential feature of the axion theory is that, when the 
axion field 9{r,t) is constant, it plays no role in electro- 
dynamics; this follows because couples to a total deriva- 
tive, e a ^ 5 F a/3 F l5 = 2e a ^ s d a {ApF 7S ), and so does not 
modify the equations of motion. However, the presence 
of the axion field can have profound consequences at sur- 
faces and interfaces, where gradients in 9(r) appear. 

A second essential feature is that electrodynamics is 
invariant under 6 — * 6 + 2~k [J]. In order to recon- 
cile this peculiar fact with the phenomenology of the 



magnetoelectric effect, observe that the axion coupling 
can alternatively be described in terms of a surface Hall 
conductivity oh whose value 9e 2 /2irh is determined by 
bulk properties, but only modulo the quantum e 2 /h. 
More generally, at an interface between two samples, 
oh = (#i — 62 + 2nr)e 2 /2%h, where the integer r de- 
pends on the details of the interface. Recall that, in 
general, a 2D gapped crystal has an integer TKNN in- 
variant C in terms of which the its Hall conductivity 
is cth = Ce 2 /h 5]. The "modulo e 2 //i", or integer r, 
discussed above corresponds to modifying the surface or 
interface by adsorbing a surface layer of nonzero C . 

When time-reversal (T) invariance is present, the 
TKNN invariants vanish, but other invariants arise that 
have been the focus of much recent work. In 2D there 
is a Z2 invariant distinguishing "ordinary" from "Z2- 
odd" insulators, with "quantum spin Hall" states 0, Q 
providing examples of the latter. In 3D there is a similar 
invariant 

S [13,113 that 

can be computed cither from the 
2D invariant on certain planes [9( or from an index in- 



volving the eight T-invariant momenta 11( . If this is odd, 
the material is a "strong topological insulator" (STI) . In 
the context of the OMP, note that T maps 9 — > —9; the 
ambiguity of 9 modulo 2-7T then implies that T invariance 
is consistent with either 9 = or 9 = tt, with the latter 
corresponding to the STI 12j . Note that if T-invariance 



extends to the surfaces, these become metallic by virtue 
of topologically protected edge states, as observed exper- 



imentally for the Bio.gSbo.i system 13J. If the surface is 
gapped by a T-breaking perturbation, then op = e 2 /2h 
modulo e 2 /h at the surface of a STI 
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In the noninteracting Berry-phase expression 

for 9 has been given in terms of the bulk bandstructure by 
Qi, Hughes, and Zhang [l2| by integrating out electrons 
in one higher dimension. Defining the Berry connection 



where 



is the cell-periodic Bloch 
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function of occupied band v and dj — d/dkj, they obtain the y direction in a supercell of length l y , they obtain 
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where the trace is over occupied bands. Note that 
wavevector-dependent unitary transformations ("gauge 
transformations") on the set of occupied wave fuctions 
cannot affect bulk physical properties. 

In the present letter, we first provide an alternate 
derivation of Eq. (0} for the OMP. Our derivation clarifies 
that 9 is a polarizability and in fact describes a contri- 
bution to magnetoelectric polarizability from extended 
orbitals. The derivation follows from an extension [15| of 
the Berry-phase theory of polarization to the case of 
slow spatial variations of the Hamiltonian. (Indeed, the 
OMP angle 9 is a bulk prop erty in exactly the same sense 
as electric polarization lall7[.) We find that the OMP 



can be generalized to the interacting case and calcu- 
lated from the many-particle wavefunction, even though 
Eq. is not valid; this reflects a subtle difference be- 
tween OMP and polarization. Explicit numerical calcu- 
lations on model crystals are presented to validate the 
theory, establish the equivalence of Eq. ((4|) to the prior 
definition, and illustrate how a non-zero 9 corresponds 
to a "fractional" quantum Hall effect at the surface of a 
magnetoelectric or topological insulator 0, [I^, EH • 

From Eq. (pj it is evident that the OMP can be viewed 
in several ways, (i) It describes the electric polarization 
arising from the application of a small magnetic field, 
(ii) It describes the orbital magnetization arising from 
the application of a small electric field, (iii) It also gives 
the (dissipationless) surface Hall conductivity an at the 
surface of the crystal, provided that the surface is insu- 
lating. Note that (iii) follows from (ii): for a surface with 
unit normal n and electric field E, the resulting surface 
current K = M x n is proportional to E x n. There is 
an elegant analogy here to the case of electric polariza- 
tion, where the surface charge of an insulating surface is 
determined, modulo the quantum e/ S, by the bulk band- 
structure alone (S is the surface cell area). 

The above discussion suggests two approaches to de- 
riving a bulk formula for the OMP 9. One is to follow 
(ii) and compute the orbital magnetization [l9| in 
an applied electrical field. We focus here on (i) instead, 
working via dP/dB. The modern theory of polarization 
starts from the polarization current jp = dP/dt under 
slow deformation of the Bloch Hamiltonian, and contains, 
to first order in d/dt, one power of the Berry curvature 
defined below [16[. Using semiclassical wavepacket dy- 
namics, Xiao et al. [l5| have shown how to compute the 
polarization current to second order and to incorporate 
slow spatial variations in the electronic Hamiltonian. For 
the case of an orthorhombic 3D crystal with M occupied 
bands in which the slow spatial variation occurs along 



for the change in the supercell-averaged polarization aris- 
ing from adiabatic currents that are inhomogeneously 
induced as a global parameter A evolves from to 1. 
Here indices ijkl run over (k x ,k y ,y,\), Tij = diAj — 
djAi — i[Ai,Aj] is the Berry curvature tensor (A^" = 
i{u^\d\\u v )) 1 and the trace and commutator refer to band 
indices. 

Because T is gauge-covariant, the integrand in Eq. ([5]) 
is explicitly gauge-invariant ; it is the non-Abelian second 
Chern class [20j , so that Eq. ([5]) is path- invariant modulo 
a quantum e/a z l y , where a z is the lattice constant in the 
z direction. Moreover, the A integral can be performed to 
obtain an expression in terms of the non-Abelian Chern- 
Simons 3- form 20]. Thus, 

(Pi in) ) =ef ^ f 1 " eij k Tr[AidjAk - 1 AA A k ] 
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where ijk now run only over (k x ,k y ,y). Here the inte- 
grand is not gauge- invariant, but the integral is gauge- 
invariant modulo the quantum e/a z l y . 

We apply this result to study the polarization 
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induced by a magnetic field described by the inhomoge- 
neous vector potential A = Byz with B — h/ea z l y , i.e., 
a B-field along x with one flux quantum threading the 
supercell. This has the effect of taking k z — > k z + eBy/h, 
and this is the only y-dependence in the Hamiltonian, so 
that \d y u) — {Be/h)\dk z u) and where ijk now run over 
(k x ,k y , k z ). Using Eqs. (|l|2j) we arrive directly at Eq. ((4]). 

There is an important geometrical relationship in this 
(nonintcracting) derivation that applies equally well to 
the many-body case and gives a bulk interpretation of 
the 2tt ambiguity in 9, whose surface interpretation was 
in terms of allowed surface IQHE layers. Polarization in a 
crystal is defined modulo the "quantum of polarization" 
flq which, for the flux-threaded supercell of Eq. ([7]), is 
AP X — e/a z l y . Since the magnetic field is B x = h/ea z l y , 
it follows that A(P X /B X ) — e 2 /h. Hence the unit-cell- 
independent ambiguity of dP/dB results from the rela- 
tionship in a finite periodic system between the unit-cell- 
dependent polarization quantum and the quantization of 
applied flux, and this relationship remains valid in the 
many-body case. 

Before studying the OMP in a specific model, we dis- 
cuss its symmetry properties and how to obtain it when 
Bloch states are unavailable, as in the many-particle case. 
Clearly the combination E ■ B in Eq. is odd under T 
and under inversion P (although it is even under the 
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combination PT). It is also odd under any improper ro- 
tation, such as a simple mirror reflection. This implies 
that 9 — —9 if the crystal has any of the above symme- 
tries. This would force an aperiodic coupling to vanish, 
but since 9 is only well-defined modulo 2tt, it actually 
only forces 9 — or tt. Thus, one can obtain an insulator 
with quantized 9 = 7r not only for T-invariant systems 
(regardless of whether they obey inversion symmetry), 
but also for inversion- and mirror-symmetric crystals re- 
gardless of T symmetry [3]. When none of these sym- 
metries are present, one generically has a non-zero (and 
non-7r) value of 9, but still retaining the simple scalar 
form of Eq. {3}. 

In an interacting system, the OMP should be obtained 
from the many-particle wavefunction. However, modify- 
ing Eq. ([J]) to the Abelian Chern-Simons integral over the 
many-body wavefunction fails 25], in important contrast 
to the case of the polarizationfthe integral of .4) , where 
such a generalization works [12]. Instead, the OMP can 
be found using the change in the many-body polarization 
due to an applied magnetic field to compute dP/dB, i.e., 
the many-body version of the supercell dP/dB calcula- 
tion. This fact is important beyond computing 9 with 
interactions, as it defines the topological insulator phase 
in the many-body case more simply than before 2l| . Like 
the IQHE, the topological insulator is defined via a re- 
sponse function (dP/dB) to a perturbation that, in the 
limit of a large system with periodic boundary condi- 
tions, is locally weak and hence does not close the insu- 
lating gap. In the IQHE, this response function is to a 
boundary phase (i.e., a flux that does not pass through 
the 2D system), while for the topological insulator, the 
defining response is to a magnetic flux through the 3D 
system. 

In the remainder of this Letter, we demonstrate the 
above theory via numerical calculations on a tight- 
binding Hamiltonian that generates non-zero values of 9, 
then discuss experimental measurements of 9. We start 
with the model of Fu, Kane, and Mele 11] for a 3D topo- 
logical insulator on the diamond lattice, 
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In the first term, the nearest-neighbor hopping amplitude 
depends on the bond direction; we take tjj = 3t + S for 
direction [111] (in the conventional fee unit cell of lin- 
ear size a) and tjj = t for the other three bonds. The 
second term describes spin-dependent hopping between 
pairs of second neighbors {{ij)), where djj and df- are 
the connecting first-neighbor legs and <x are the Pauli 
spin matrices. With \S\ < It and Xso sufficiently large, 
this model has a direct band gap of 2|<5|. 

To break T we add a staggered Zeeman field with 
opposite signs on the two fee sublattices A and B, 

h ' (j2ieA 4 trc i ~ J2ieB c \ <TC i) ■ We take \ h \ = "isin^ 



and choose h in the [111] direction; setting 6 = mcos/3 
and varying the single parameter (3 keeps the gap con- 
stant and interpolates smoothly between the ordinary 
((3 = 0) and the topological ((3 = n) insulator. 

We have calculated the OMP angle 9 using four dif- 
ferent methods with excellent agreement (Fig.[T]). First, 
we obtain 9 from Eq. ((4]); this requires a smooth gauge 
for A, which can be found using now-standard Wannier- 
based methods [51]. Results are shown for — ir/4 and 
(3 = tt/2 (filled squares). 
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FIG. 1: The magnetoelectric polarizability 9 (in units of 
e 2 /2ivh). The filled squares are computed by the Chern- 
Simons form, Eq. ((4]). The open squares are dP/dB from 
Eq. ((9]l . The points are obtained by layer-resolved oh calcu- 
lations using Eq. (|12[) . The curve is obtained from Eq. (|10[) . 



Next, we have calculated the polarization [if 
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resulting from a single magnetic flux quantum in a large 
supercell. Varying the supercell size (and thereby B) al- 
lows us to approximate dPj dB, yielding the open squares 
in Fig. Q] The points in Fig. [T] are from the surface Hall 
response in a slab geometry, described below. Finally, to 
obtain the curve in Fig. [TJ we also computed 0(0) from 
the second Chern expression U HI 



d0 / d 3 fce^Tr[^(k,/3')^(k,/?')] (10) 
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(derived above as Eq. |JSJ|). Clearly, the various ap- 
proaches are numerically equivalent. 

We now discuss the surface Hall conductivity, whose 
fractional part in units of e 2 /h is just 9/2ir A]. Consider 
a material with coupling 9 in a slab geometry that is finite 
in the z direction and surrounded by 9 = vacuum. The 
simplest interfaces will then lead to <jh = 0e 2 /(2irh) at 
the top surface and —9e 2 /(2Trh) at the bottom surface, 
for a total a xy of zero. More generally, arbitrary surface 
quantum Hall layers change the total integer quantum 
Hall state, but not the fractional parts at each surface. 

The spatial contributions to the Hall conductance in 
the slab geometry can be resolved as follows. The unit 
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FIG. 2: (Color online) The layer-resolved Hall conductivity 
(in units of e 2 /h) at f3 — ir in a slab of twenty layers, with 
m = t/2 and Ago = t/i, terminated in (111) planes. 



Landau levels may make it difficult to isolate the half- 
integer part of surface cth- In the presence of broken 
discrete symmetries, as in antiferromagnets or multifer- 
roics, the surface gap exists naturally and experiments 
are easier. For example, the theoretical methods of this 
paper could be used to compute the orbital part of the 
recently measured 9 in Cr 2 03 

The authors acknowledge useful discussions with A. Se- 
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0804413 (JEM), and NSF DMR-0549198 (DV). 



cell is a supercell containing some number N of original 
unit cells in the z direction, with translational invariancc 
remaining in the x and y directions. The TKNN integer 
for the entire slab is [B|, |23| 
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Here i and j take the values k x and k y and V — 
y] t , is the projection operator onto the occupied 

subspace (y runs over occupied bands). To find how dif- 
ferent i layers contribute to C, define a projection V n 
onto layer n within the supercell, and compute 



C{n) = -k 
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The results, presented in Fig. [2j confirm that the surface 
layers have half-integer Hall conductance when [3 = it in 
((HI) and that the sign on each surface is switched by lo- 
cal T-breaking perturbations (in this example, a uniform 
Zeeman coupling in the surface layer). 

To gain some insight into the microscopic origin of 9 in 
the noninteracting case, using Eq. we have calculated 
9 for a Hamiltonian that breaks PT (as well as P and T) 
by adding a weak, uniform (i.e., not staggered) Zeeman 
coupling. For some values of (3 this lifts all degenera- 
cies, enabling us to isolate the single-band and interband 
contributions to 9 and to verify that, because interband 
contributions are nonzero in general, 9 is a property of 
the whole occupied spectrum (unlike polarization, which 
is a sum of individual band contributions). A single filled 
band can have nonzero 9 only if there are more than two 
bands in total 



24] 



Experimental detection of 9 is more difficult for a topo- 
logical insulator than for a generic magnetoelectric insu- 
lator because some T-breaking perturbation is needed 
to gap the surface state. Furthermore, a large surface 
density of states, as in Bio.gSbo.i, may complicate the 
measurement: while even a weak magnetic field will in 
principle lead to a gap and half-integer quantum Hall ef- 
fect at each surface, the large number of filled surface 



[1] N. A. Spaldin and M. Fiebig, Science 309, 391 (2005). 
[2] M. Fiebig, J. Phys. D - Applied Phys. 38, R123 (2005). 
[3] F. W. Hehl, Y. N. Obukhov, J.-P. Rivera, and H. Schmid, 

Physics Letters A 372, 1141 (2008), 0708.2069. 
[4] F. Wilczek, Phys. Rev. Lett. 58, 1799 (1987). 
[5] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and 

M. den Nijs, Phys. Rev. Lett. 49, 405 (1982). 
[6] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 226801 

(2005). 

[7] C. L. Kane and E. J. Mele, Phys. Rev. Lett. 95, 146802 
(2005). 

[8] B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science 

314, 1757 (2006). 
[9] J. E. Moore and L. Balents, Phys. Rev. B 75, 121306(R) 

(2007). 

[10] R. Roy, cond-mat/0607531. 

[11] L. Fu, C. L. Kane, and E. J. Mele, Phys. Rev. Lett. 98, 
106803 (2007). 

[12] X. Qi, T. L. Hughes, and S.-C. Zhang, Phys. Rev. B 78, 
195424 (2008). 

[13] D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. 

Cava, and M. Z. Hasan, Nature 452, 970 (2008). 
[14] L. Fu and C. L. Kane, Phys. Rev. B 76, 045302 (2007). 
[15] D. Xiao, J. Shi, D. P. Clougherty, and Q. Niu, Phys. Rev. 

Lett. 102, 087602 (2009). 
[16] R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 

1651 (1993). 

[17] G. Ortiz and R. M. Martin, Phys. Rev. B 49, 14202 
(1994). 

[18] T. Thonhauser, D. Ceresoli, D. Vanderbilt, and R. Resta, 

Phys. Rev. Lett. 95, 137205 (2005). 
[19] D. Xiao, J. Shi, and Q. Niu, Phys. Rev. Lett. 95, 137204 

(2005). 

[20] M. Nakahara, Geometry, Topology and Physics (Institute 

of Physics Publishing (Bristol), 1998). 
[21] S.-S. Lee and S. Ryu, Physical Review Letters 100, 

186807 (2008), 0708.1639. 
[22] N. Marzari and D. Vanderbilt, Phys. Rev. B 56, 12847 

(1997). 

[23] J. E. Avron, R. Seiler, and B. Simon, Phys. Rev. Lett. 
51, 51 (1983). 

[24] J. E. Moore, Y. Ran, and X.-G. Wen, Physical Review 
Letters 101, 186805 (2008), 0804.4527. 

[25] Mathematically, obtaining a nontrivial second Chern or 
Chern-Simons integral depends on having a degenerate 
set of bands somewhere in parameter space, which is not 
the case for a gapped many-body wavefunction. 



